# ------------------------------------------------------------------------------
# Replication Materials
# 
# title: Eliciting Beliefs as Distributions in Online Surveys
# journal: Political Analysis
# authors: Lucas Leemann, Richard Traunmüller, and Lukas Stoetzer
# date: August 2020
# ------------------------------------------------------------------------------


  library(tidyverse)
  library(xtable)
  options(xtable.comment = FALSE)
  source("fun/fun_electbeleifs_MLE.R")
  source("fun/fun_utilities.R")

  
  

  
  
  
  source("00_DataCleaner_addon.R") # Richard Code to clean Data
  dat_sym <- filter(data, screen1==1) # Filter with screen questions


# Estimation  ============

  res <- list() # Result List


## Quantile Question  ============

  # Estimate Model for Quantile Question (with correction?)
  res[["quant"]] <- dat_sym %>%
    creat_dat_Y(method="Quantile") %>%
    est_eB_beta(true.val = c(1,1))
  
  res[["quant.dem"]] <- dat_sym %>%
    filter(PID=="Democrat") %>%
    creat_dat_Y(method="Quantile") %>%
    est_eB_beta(true.val = c(1,1))

  res[["quant.rep"]] <- dat_sym %>%
    filter(PID=="Republican") %>%
    creat_dat_Y(method="Quantile") %>%
    est_eB_beta(true.val = c(1,1))
  

      
## Intervall Question (Wide) ============

  # Estimate Model for WIDE Interval
  res[["intwi"]]  <- dat_sym %>%
    creat_dat_Y("Interval") %>%
    est_eB_beta_int(cut.offs = c(0.40,0.60), true.val = c(1,1))
  
  res[["intwi.dem"]]  <- dat_sym %>%
    filter(PID=="Democrat") %>%
    creat_dat_Y("Interval") %>%
    est_eB_beta_int(cut.offs = c(0.40,0.60), true.val = c(1,1))
  
  res[["intwi.rep"]]  <- dat_sym %>%
    filter(PID=="Republican") %>%
    creat_dat_Y("Interval") %>%
    est_eB_beta_int(cut.offs = c(0.40,0.60), true.val = c(1,1))

  
## Intervall Question (Narrow) ============

  # Estimate Model for NARROW interval
  res[["intna"]]  <- dat_sym %>%
    creat_dat_Y("Interval2") %>%
    est_eB_beta_int(cut.offs = c(0.45,0.55), true.val = c(1,1))

  res[["intna.dem"]]  <- dat_sym %>%
    filter(PID=="Democrat") %>%
    creat_dat_Y("Interval2") %>%
    est_eB_beta_int(cut.offs = c(0.45,0.55), true.val = c(1,1))
  
  res[["intna.rep"]]  <- dat_sym %>%
    filter(PID=="Republican") %>%
    creat_dat_Y("Interval2") %>%
    est_eB_beta_int(cut.offs = c(0.45,0.55), true.val = c(1,1))
  
  
## Manski Question  ============

  # Estimate Model for Manski
  res[["manski"]]  <- dat_sym %>%
    creat_dat_Y("Manski") %>%
    est_eB_beta_manski(true.val = c(1,1))

  res[["manski.dem"]]  <- dat_sym %>%
    filter(PID=="Democrat") %>%
    creat_dat_Y("Manski") %>%
    est_eB_beta_manski(true.val = c(1,1))
  
  res[["manski.rep"]]  <- dat_sym %>%
    filter(PID=="Republican") %>%
    creat_dat_Y("Manski") %>%
    est_eB_beta_manski(true.val = c(1,1))
  

  
## Bins and Balls Question  ============

  # Estimate Model for Bins & Balls
  res[["binba"]]  <- dat_sym %>%
    creat_dat_Y("BinsBalls") %>%
    est_eB_beta_binsballs(true.val = c(1,1))
  
  res[["binba.dem"]]  <- dat_sym %>%
    filter(PID=="Democrat") %>%
    creat_dat_Y("BinsBalls") %>%
    est_eB_beta_binsballs(true.val = c(1,1))
  
  res[["binba.rep"]]  <- dat_sym %>%
    filter(PID=="Republican") %>%
    creat_dat_Y("BinsBalls") %>%
    est_eB_beta_binsballs(true.val = c(1,1))

#+ Figure, echo=FALSE, include=TRUE, fig.cap = "Comparsion of average elecited beliefs and true data distribution"
## Figure Average Estimate  ============

  # Get estimates
  par_est <- data.frame(t(sapply(sapply(res, "[[", "par"),
                                 function(res) res[names(res) %in% c("alpha","beta")])))


  # Expand Data.frame
  expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))
  df_est <- expand.grid.df(data.frame("x"=seq(0.2,0.8,0.01)), par_est)
  df_est$format <- factor(kronecker(c(1:15),rep(1,61)), labels = c("Quantile", "DEM: Quantile","REP: Quantile",
                                                                  "Interval (Wide)", "DEM: Interval (Wide)", "REP: Interval (Wide)",
                                                                  "Interval (Narrow)", "DEM: Interval (Narrow)", "REP: Interval (Narrow)",
                                                                  "Manski","DEM: Manski","REP: Manski",
                                                                  "Bins and Balls","DEM: Bins and Balls","REP: Bins and Balls"
                                                                  ))
  df_est$Sample <- factor(kronecker(c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3), rep(1,61)), labels=c("All","Democrats","Republicans"))
  df_est$Sample1 <- factor(abs(kronecker(c(1,0,0,1,0,0,1,0,0,1,0,0,1,0,0), rep(1,61))-1), labels=c("All", "by Party ID"))
  df_est$plot.num <- factor(kronecker(c(2,2,2,4,4,4,5,5,5,1,1,1,3,3,3),rep(1,61)), 
                                      labels=c("Manski","Quantile", "Bins and Balls","Interval (Wide)", "Interval (Narrow)"))
  
    
  # Prepare Data for Plot
  df_est  <- mutate(df_est, db = dbeta(x,alpha,beta))


  # Figure 
  ggplot(df_est, aes(x=x, y=db, group=Sample, size = Sample, color=Sample)) +
    geom_line() +
    facet_grid(plot.num~Sample1) +
    theme_minimal() +
    scale_color_manual(values = c("black","blue","red")) +
    scale_size_manual(values = c(1.5,1.5,1.5)) +
    xlab("Vote Share") + ylab("Density") +
    theme(axis.title.y=element_blank(),
          text = element_text(size=20),
          strip.text.y = element_text(angle = 360),
          legend.position = "bottom",
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          panel.background = element_rect(fill = "grey96",
                                          colour = "grey96",
                                          size = 0, linetype = "solid"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank() 
    )

  ggsave("out/fig/Fig_5.pdf", width=12, heigh=10)
  
# add for table
  par_est1  <- par_est %>%
    mutate(mean = alpha/(alpha + beta)) %>%
    mutate(variance = alpha*beta/((alpha+beta+1)*(alpha+beta)^2))
  
  par_est2 <- par_est1
  par_est2[4:6,] <-  par_est1[1:3,]
  par_est2[10:12,] <-  par_est1[4:6,]
  par_est2[13:15,] <-  par_est1[7:9,]
  par_est2[1:3,] <-  par_est1[10:12,]
  par_est2[7:9,] <-  par_est1[13:15,]
  
  
  
  
  print(xtable(par_est2, caption="Parameter Estimates & Moments of Belief Distribution"), floating = FALSE,booktabs = T, include.rownames=FALSE,
        file=paste("out/fig/Tab_A10.tex",sep=""))
  
  

